Glassy transition and met ast ability in four-spin Ising 

model 

Adam Lipowski 1),2) and Des Johnston 1 ) 

1} Department of Mathematics, Heriot-Watt University, 
EH14 4AS Edinburgh, United Kingdom 
^Department of Physics, A. Mickiewicz University, 
61-614 Poznan, Poland 

Abstract 

Using Monte Carlo simulations we show that the three-dimensional Ising model with four- 
spin (plaquette) interactions has some characteristic glassy features. The model dynamically 
generates diverging energy barriers, which give rise to slow dynamics at low temperature. 
Moreover, in a certain temperature range the model possesses a metastable (supercooled liq- 
uid) phase, which is presumably supported by certain entropy barriers. Although extremely 
strong, metastability in our model is only a finite-size effect and sufficiently large droplets 
of stable phase divert evolution of the system toward the stable phase. Thus, the glassy 
transitions in this model is a dynamic transition, preceded by a pronounced peak in the 
specific heat. 

§ 1. Introduction 

Recently, glassy systems have been intensively studied from both the theoretical and exper- 
imental points of view |2[. However, due to the very complicated nature of glasses their 
understanding is still far from being complete. Although one can construct various off-lattice 
models of glasses, these models are in general very difficult to study and additional insight would 
be desirable Q . 

One possible direction is to study lattice models of glasses in the hope that, albeit unreal- 
istic in some respects, these models describe at least some aspects of a glassy transition. The 
advantage of using lattice models is that usually they are much easier to study. Indeed, one can 
show analytically, and in some cases even exactly, that certain lattice models do undergo glassy 
transitions 0j. However, these analytically tractable models are usually of infinite dimension 
or contain interactions of infinite range and it is not certain to what extent these models are 
applicable to real systems. 

An interesting model exhibiting some glassy features, which is both three-dimensional and 
contains short-range interactions, was proposed some years ago by Shore et al. |6|. They 
studied the dynamics of a three-dimensional ferromagnetic Ising model with antiferromagnetic 
next-nearest-neighbours interactions (the SS model). They showed that the low-temperature 
coarsening of a random quench asymptotically becomes very slow and the characteristic length 
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scale increases logarithmically in time, which is one of the characteristic features of glasses. 
However, Shore et al. argued that their model is not yet a satisfactory model of glasses. Namely 
they show that in their model the transition into the slow-dynamics regime is induced by a corner- 
rounding transition. At this transition the energy barriers should vanish, which, according to 
their calculations, implies that the zero-temperature characteristic length scale increases too fast 
as a function of inverse cooling rate ||. Shore et al. argued that in more realistic models of 
glasses the energy barriers should exist even at the glassy transition. 

In our opinion, the SS model lacks yet another characteristic feature of real systems, namely 
the supercooled liquid phase. In the SS model the high-temperature phase (liquid) cooled below 
the critical point T c but above the corner-rounding transition relatively quickly evolves toward 
the low-temperature phase (crystal). One would expect that when brought to this temperature 
range, liquid should not immediately crystallize (or rather polycrystallize) but, at least for some 
time, it should remain in a metastable state of supercooled liquid. 

Recently, it has been shown that the three-dimensional Ising model with four-spin (plaquette) 
interactions has some similarities to the SS model |7]]. Namely, it was shown that in this model 
certain configurations are very long-lived due to large energy barriers. Configurations of this kind 
should develop spontaneously during the coarsening and thus should considerably slow-down this 
process. Indeed, it was observed that the high-temperature quench becomes trapped in a 
glassy phase. 

In the present paper we examine this model further. First, we show that in the low- 
temperature regime the characteristic length l(t) increases very slowly in time. We argue that 
this increase is likely to be logarithmic, namely l(t) ~ Int. Such a slow increase of the char- 
acteristic length is a typical feature of glasses. Although the origin of the very slow dynamics 
in our model is basically the same as in the SS model, the nature of the glassy transition is 
quite different. Firstly, above the glassy transition the model does not enter the fast-dynamics 
regime (as the SS model does) but remains trapped in the supercooled liquid phase. In addition, 
our estimation of a certain characteristic time shows that energy barriers exist even above the 
glassy transition. It strongly suggests that the glassy transition in our model is not induced by 
a corner-rounding transition and that energy barriers persist even above the glassy transition. 

There are other properties of the model which are very interesting. Our simulations show 
that in this model in a certain temperature range, depending on the initial conditions, the model 
might be either in the liquid or the crystal phase. Using the thermodynamic integration method 
we calculated the free energies of both phases and, as expected, the crossing point (i.e., the 
first-order phase transition) appears approximately in the middle of this temperature range. 
Such metastable effects (hysteresis) frequently accompany first-order phase transitions || |9j. A 
distinctive feature of metastability of our model is that this is an extremely strong effect. We 
show that even for temperatures close to the limits of hysteresis only a large droplet of the stable 
phase injected into the metastable phase can divert the evolution of the system into the stable 
phase. Spontaneous nucleation of such large droplets is an extremely improbable event, and is 
well beyond computational capacities of modern computers. 

In addition, our simulations show that upon cooling the glassy transition is accompanied by 
a pronounced peak in the specific heat. This (pseudo-)critical behaviour gives rise to certain 
slow modes in the liquid phase of our model. Such slow modes for glassy systems are predicted 
by the Mode-Coupling Theory j!J and are experimentally verified as so-called a-oscillations [ [lC[| . 

Our paper is organised as follows. In section 2 we introduce the model, review its basic 
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properties and study the evolution of random quench. In section 3 we calculate the free energy 
and the specific heat and present the time evolution of internal energy. In section 4 we discuss 
metastability properties of our model. Section 5 contains our conclusions. 



§ 2. Basic properties and domain coarsening 

The model is defined by the following Hamiltonian 

H = -J2S i S j S k S h (2.1) 

where Si = ±1 and summation is performed over all elementary plaquettes k, I) of a simple 
cubic lattice of linear size L with spins placed at the sites of the lattice. Model ( |2.1| ) is a special 
case of the so-called gonihedric models, which have recently been studied in the context of the 



lattice field theory [11, |12|| . Closely related models have been studied in the context of random 
surface models and a very rich phase structure mapped out |i"4|| . 

The basic properties of model (^T|) which are already known can be briefly described as 
follows: The ground state is strongly degenerate with the degeneracy ~ 2 3L (the ground state 
entropy per site is thus zero). Monte Carlo simulations show [|l^, [?]] that upon heating of 
an arbitrary ground state (crystal) configuration, the model undergoes a discontinuous phase 
transition into the disordered phase (liquid) around T ~ 3.9; the temperature scale is set as 
in 0] with Boltzmann constant put to unity. An important feature of model (^^) is the shape 
dependence of excitations J?]: it is not only the size of the excitation which determines its energy 
(as in the ordinary two-spin Ising model) but also its shape (see Fig. [l]). This property, which 
holds also for the SS model, gives rise to energy barriers which are in turn responsible for the 
very slow dynamics of both models. 

To study the evolution of the random configuration quenched to low temperature, we mea- 
sured the energy excess SE(t) = E(t)—Eo over the ground state energy Eq = —3. One expects || 
that the inverse of this quantity sets the characteristic length scale l(t) of the system, which 
roughly corresponds to the average size of domains. Moreover, there is convincing evidence [|ll| 
that in many systems with nonconservative dynamics and a scalar order parameter (i.e., condi- 
tions which are satisfied in our approach) l(t) increases asymptotically in time as l{t) ~ t n and 
n = 1/2. However, in some systems l{t) is known to increase much more slowly in time, namely 
logarithmically l(t) ~ log(i). These exceptional systems include some random (at the level of the 
Hamiltonian) systems []i~6[ |i~7| , and the SS model for temperatures below the corner-rounding 
transition flg|| . It is the energy barriers developing in these systems during the coarsening which 
cause such a slow increase of l(t). 

We performed standard p3] Monte Carlo simulations using the Metropolis algorithm with a 
random sequential update. Unless stated otherwise, periodic boundary are imposed. The log-log 



plot of l/6E(t) as a function of time for model (2.1) is shown in Fig. |2|. The presented results 



are obtained for L = 40 but very similar behaviour was observed for L = 30. From Fig. Q it is 
clear that for T = 1.5 and 2.8 the asymptotic slopes of the curves are much smaller than 1/2 
and there is a tendency for these curves to bend downwards. Taking into account the absence of 
models with n considerably smaller than 1/2 and the existence of energy barriers in model ( |2.1| ) 
of basically the same nature as in the SS model, suggests that for the examined temperatures 
the characteristic length asymptotically increases logarithmically in time. We cannot exclude, 
however, that in this temperature regime the increase of l{t) is even more exotic, with neither 
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logarithmic nor power-law increase. It would appear that such a slow increase of lit) takes place 
even for T = 3.3 and 3.4, but the behaviour of l{t) for these temperatures is obscured by the 
metastability effects, since before collapsing into the glassy phase the system remains in the 
liquid state for some time. 

The difference between our model and the SS model becomes clear when we approach the 
glassy transition (which we roughly estimate to take place at T = T g ~ 3.4) by increasing the 
temperature. In the SS model for temperatures below the critical point but above the corner- 
rounding transition thermal fluctuations roughen corners of domains and energy barriers are 
no longer relevant. Consequently, the "ordinary " dynamics with n = 1/2 is restored and the 
system rapidly evolves toward the low-temperature phase. On the contrary, in model ( |2.1| ) for 
3.4 < T < 3.9 the random quench does not even evolve toward the low-temperature phase but 
remains disordered Since the low-temperature phase does exist at these temperatures (as 
we have already mentioned, the transition from the low-temperature phase to the disordered 
phase takes place around T ~ 3.9), there should be some barriers which prevent the liquid from 
collapsing. In our opinion, these barriers are of entropic origin and they are probably related to 
the strong degeneracy of the ground state, which would explain why the behaviour of our model 
is different from the SS model [ 20 [ . This phenomenon is discussed further in section 4. 



§ 3. Free energy and specific heat 

The free energy encodes all the important thermodynamic information about a system. In 
equilibrium statistical mechanics this quantity is defined as 

f = -T/Nln[J2e- H / T ], (3.1) 

where N is the number of particles (lattice sites) and H is the Hamiltonian of the system with 
the Boltzmann constant put to unity. However, calculation of the above defined free energy using 
Monte Carlo simulations is not straightforward and requires thermodynamic integration [pJ. To 
calculate the free energy of our model in the liquid and crystal phases we used the following 
equations: 

f T c [ X I T 1 

/cryst = u-T —dT, /iiq = -Ts(oo) + T / ud(— ). (3.2) 
Jo J Jo 1 

In the above equation c and u denote the specific heat and internal energy, respectively 
(calculated using standard formulae |]l9|), and s(oo) = ln(2) is the entropy per site at infinite 
temperature. The results of our calculations are shown in Fig. [3|. We checked the stability of our 
results with respect to the integration step AT (or A(l/T)), the system size L and the number 
of Monte Carlo steps at each temperature. 

In Fig. U one can see that the free energies of the crystal and of the liquid cross around T = 3.6 
and we expect that this is the temperature of the first-order phase transition. The estimation 
of the transition temperature is in a good agreement with calculations done using the Cluster 
Variational Method [plj| . However, due to the strong metastability, which is discussed in more 
detail in the next section, the transition at this point is very difficult to observe. Indeed, if 
we heat a crystal sample, the transition occurs around T = 3.9, while cooling a liquid sample 
results in the glassy transition around T = 3.4. 

To overcome the metastability and confirm that a first-order transition does take place around 
T = 3.6, we simulated the system with a nonuniform initial configuration. Namely, we prepared 



4 



the system with one half of it in the crystal phase (e.g., all spins 'up') and the second half in the 
liquid phase (all spins random); see Fig. |4j. Such a choice facilitates the evolution toward the 
stable phase (i.e., the one with the lowest free energy), because both phases are present in the 
initial configuration and the system does not have to nucleate the stable phase. Results of such 
simulations are shown in Fig. |5[ One can clearly see that the evolution of the system depends 
on whether the temperature is above or below T = 3.6. For T < (>)3.6 the crystal (liquid) 
phase gradually expands until the stable phase invades the whole system. The identification of 
the final state is obtained from a comparison of its internal energy with simulations which use a 
uniform initial state and also from the visual inspection of snapshot Monte Carlo configurations. 

To provide additional information about model (2.1), we measured the variance of the in- 
ternal energy and calculated the specific heat |19[| . Our data is shown in Fig. ^. When we start 
our simulations from the ground-state configuration (heating), the behaviour of the specific heat 
confirms a transition around T = 3.9, in agreement with earlier simulations fll2| , [7j. However, 
under cooling this peak is shifted toward much lower temperature and it seems to coincide with 
our estimation of T g . Under cooling we do not observe any singularity in the specific heat until 
T = T g , which indicates that during cooling and for T > T g the system remains in the liquid 
phase. Let us notice that the locations of both peaks are almost unchanged after increasing the 
system size by almost a factor of two. Moreover, for T < T g the specific heat is slightly larger 
upon cooling than upon heating. This is in agreement with the fact that in this temperature 
range the model has slow dynamics and cannot reach (within numerically accessible computing 
time) the crystal phase, for which the specific heat is very small. 

Let us note that the pronounced peak in the specific heat presumably indicates critical or 
pseudo-critical behaviour in model (|2.1[ ). If so, we might expect that certain relaxation times 
might substantially increase or even diverge, which should be observed as some slow modes in 
the system. In Fig. |7|we present the time dependence of the energy E(t). Although an explicit 
calculation of the Fourier transform of this quantity would be desirable, it is rather clear that 
for T = 3.42, which is very close to T g , in addition to fast fluctuations the system exhibits 
slow fluctuations with a time scale ~ 150. For higher temperatures (T = 3.6) the time scale 
of slow fluctuations decreases and eventually (T = 3.8) becomes hardly distinguishable from 
fast oscillations. Such slow and fast modes resemble experimentally observed a— (slow) and 
j3— (fast) oscillations in real glasses [10]. It would be interesting to check to what extent the 
properties of our model agree with the Mode Coupling Theory HJ , according to which such slow 
modes are a key factor driving a glassy transition. 



§ 4. Metastability 

Results presented in the previous sections suggest that in the temperature range T g < T < T c 
the system might remain either in the liquid or the crystal phase. In the present section we 
present some additional results concerning the (meta-)stability of the liquid and crystal phases. 



4.1 Characteristic times 

We measured various characteristic times imposing different initial and boundary conditions and 
monitoring the evolution toward a final state. To check the stability of the liquid (Ti; q ), we used 
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a random initial configuration and simulated the system until the energy reached E = -2.3. [j 
To calculate Ty iq we made 100 independent runs. Our results for T = 3.5 are shown in Fig. || 
and they suggest that the escape time increases at least exponentially with the system size. 

To check the stability of the crystal, one should measure the time (r+) needed for the crystal 
to be transformed into the liquid. It would be particularly interesting to examine the size 
dependence of r + for 3.6 < T < 3.9, i.e., for temperatures where the crystal is metastable. 
We have found, however, that this quantity increases very rapidly with the system size and 
in this temperature range it is virtually impossible to increase the system size beyond L = 6. 
The stability of this phase might be also inferred from other measurements we made, in which 

we estimated the time (r_| ) needed to shrink a cubic excitation of size L. This technique 

parallels that which has already been applied to the SS model Q: the initial configuration has 
"up "spins at the boundary of the cube of size L + 2 (which are kept fixed) and "down "spins 
inside this cube. Simulations are performed until the magnetization of the interior of the cube 

decays to zero and the time needed for such a run is recorded. To calculate at a given 

temperature we made 100 independent runs. Our results for T = 3.6 (Fig. |8|) suggest that 

t _| increases approximately exponentially with L. Such an increase confirms that the glassy 

transition in model ( |2.lD is not induced by the corner-rounding transition since above the corner- 
rounding transition one expects r_| ~ L 2 || and the data in Fig. || should bend considerably 

downwards. It also confirms the stability of the crystal since it is clear that bringing the crystal 
(a homogeneous, low-energy configuration) into the liquid is a slower process than shrinking an 

excitation. Similar size dependence of rjiq and t_j was also observed for other temperatures in 

the interval 3.4 < T < 3.9. 

We also measured the characteristic times Tn q and r + outside the interval 3.4 < T < 3.9. 
Simulations were done for several values of system size L and the results were extrapolated to 



the thermodynamic limit (L — > oo) using a simple fit [22]. These extrapolated values are shown 
in Fig. |9|. One can see that outside this temperature range rii q and r + are definitely finite and 
they seem to diverge upon approaching T = 3.4 from below (rii q ) and T = 3.9 from above (t+). 

4.2 Droplets of a stable phase 

The numerical data presented in the previous subsection suggests that for 3.4 < T < 3.9 
the model has two different phases of effectively infinite life-time. Such a result would be 
in disagreement with the result that in short-range interacting systems metastability is only 
a quantitative effect |^, |9|. In this section we show, however, that data presented in Fig. || 
and Fig. ^ are misleading and in the thermodynamic limit {L — > oo) the characteristic times rii q 
for T < 3.6 and r+ for T > 3.6 should be finite. The time scale of these metastable effects is, 
however, enormously long in comparison with the length of our simulations. 

One expects that metastable phases have only a finite life-time due to droplet nucleation. 
When a sufficiently large droplet of stable phase nucleates inside a metastable phase, it diverts 
evolution of the system toward a stable phase. Since the critical (i.e., minimal) droplet size is 
finite, there is a finite probability of spontaneous nucleation of such droplets and thus a life-time 
of a metastable phase is also finite. 



To check whether such a mechanism operates in model (2.1), we monitored the evolution 



of the system with droplets introduced by hand into the initial configuration (see Fig. |j). Our 

This value is chosen rather arbitrarily, but once the system reaches this energy it does not return to the liquid 
phase. 
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results for T = 3.5 are shown in Fig. |T^. One can see that when a droplet of the crystal is 
sufficiently small (of linear size M = 18,24), the system after some transient ends up in the 
liquid phase. However, a large droplet of size M = 36 diverts evolution of the system toward 
the more stable (crystal or glassy) phase. We performed similar simulations to examine the 
(meta-)stability of the crystal phase. Setting T = 3.8 we observed that droplets of liquid of size 
M > 24 divert evolution of the crystal phase toward the stable liquid phase. 



The above results show that the droplet-nucleation mechanism is effective in model (2.1 



and metastable phases are of finite life-time. However, the important question is how long is 
this life-time. This quantity is determined by the inverse of the probability of the spontaneous 
nucleation of critical droplets. It is clear that for computationally accessible systems (L ~ 100) 
spontaneous nucleation of crystal droplets of linear size M ~ 30 is an extremely unlikely event, 
which takes place on astronomical time scales. We should emphasize that it does not mean that 
our model predicts such a life-time of metastable liquids. Since the nucleation of droplets is 
basically a local event, its probability for macroscopic systems increases merely due the system 
size (droplets might nucleate independently in many places). 

The radius of critical droplets presumably vanishes upon approaching the limits of hysteresis 
(i.e., T = 3.4 and 3.9). Thus, very close to these limits one should be able to observe finite- 
life-time effects such as the spontaneous collapse of liquid into the crystal (or maybe glassy) 
phase. 



§ 5. Conclusions 

In the present paper we have studied the three-dimensional Ising model with four-spin inter- 
actions. The Hamiltonian of this model is homogeneous, non-frustrated and contains only 
short-range (plaquette) interactions. Nevertheless, we found that this model has very interest- 
ing dynamical and thermodynamical properties. In particular, our numerical results suggest 
that the model has a very slow coarsening dynamics in its low-temperature phase. Moreover, 
due to very strong metastability, in a certain temperature range the model can remain either 
in a crystal or liquid phase depending on how it has been prepared. We have shown that 
a droplet-nucleation mechanism is effective in this model and thus that metastability in this 
model is a finite-size (but very strong) effect. The time scale of spontaneous nucleation for crit- 
ical droplets is extremely large and well beyond the timescales of our simulations. Upon cooling 
the metastable liquid collapses into a glassy phase. We have found that the glassy transition 
in our model, which is purely dynamical in nature, is preceded by a pronounced peak in the 
specific heat (pseudo-critical behaviour) and also by slow oscillations of the internal energy. 
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Figure captions 



Fig. 1: In model ( |2.1| ) the energy of an excitation (e.g., 'down' spins surrounded by 'up' spins) is 
proportional to the total length of edges of the boundary of that excitation. To remove 
a cubic excitation (a), the system is likely to proceed through configurations like that 
shown in (b) and (c). It is easy to realize that in (b) and (c) the total length of edges 
and thus the energy of such configurations is larger than that in (a). In the case (c) the 
energy increase is proportional to the linear size of the excitation. On the contrary, in a 
two-spin Ising model all configurations would have the same energy (since the area of all 
configurations is the same) and removing of excitations would proceed without climbing 
any energy barriers. 

Fig. 2: The log-log plot of of l/5E(t) as a function of time t (L = 40). The dashed line has a 
slope 1/2. 

Fig. 3: The free energy of liquid (+) and crystal (o). Calculations were done for L = 30, and 
5000 Monte Carlo steps were used at each temperature. The integration steps were AT = 
A(l/T) = 0.025. 

Fig. 4: Two-dimensional sections of initial configurations used in our calculations, (a) An initial 
configuration used in the calculations of internal energy shown in Fig. [|. (b) An initial 
configuration with a droplet of crystal phase injected into the liquid phase (see Section 
4.2). Spins in the crystal/liquid part of the system are initially set 'up'/at random. 

Fig. 5: The internal energy U as a function of time for L = 50 and an initial configuration 
composed of 50% of crystal and 50% of liquid (see Fig. |3|a). The steady-state values of 
U for T = 3.5, 3.55 and for T = 3.65, 3.7 are identical (within error bars) with internal 
energy of crystal and liquid at corresponding temperatures. 

Fig. 6: The specific heat calculated during (i) heating for L = 24 (o) and L = 40 (x) (ii) cooling 
for L = 24 (+) and L = 40 (□). For each temperature and system size we made runs of 
10 4 Monte Carlo steps plus 10 3 Monte Carlo steps for relaxation. 

Fig. 7: The time dependence of the energy E(t) in the liquid phase at various temperatures 
(L = 40). 

Fig. 8: The size dependence of the logarithm of the escape times Tii q (□) and r_| (+). Calculation 

of Tn q and t_| was done for T = 3.5 and T = 3.6, respectively. 

Fig. 9: The inverse of the characteristic times rii q (+) and r + (o) as a function of temperature. The 
plotted results are obtained by extrapolation of the finite-size data to the thermodynamic 
limit. 



Fig. 10: 



The internal energy U as a function of time for T = 3.5. The initial configuration consists 
of a droplet of crystal phase of size M injected into the liquid phase. 



